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1. Introduction 

The radial Schrodinger equation can be written as: 

y'\x) = [1(1 + l)/x^ + Vix) - k'^]y{x). (1) 

Many problems in theoretical physics and chemistry, material sciences, 
quantum mechanics and quantum chemistry, electronics etc. can be 
express via the above boundary value problem (see for example [1] - 
[4]). 

We give the definitions of some terms of (1): 

— The function W{x) = l{l + -|- V{x) is called the effective 
potential. This satisfies W{x) ^ as x ^ cxo 

— The quantity k'^ is a real number denoting the energy 
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— The quantity I is a given integer representing the angular momen- 
tum 

— F is a given function which denotes the potential. 
The boundary conditions are: 

y(0) = (2) 

and a second boundary condition, for large values of x, determined by 
physical considerations. 

The last years an extended study on the development of numerical 
methods for the solution of the Schrodinger equation has been done. 
The aim of this research is the development of fast and reliable methods 
for the solution of the Schrodinger equation and related problems (see 
for example [5] - [18], [24] - [127]). 

We can divide the numerical methods for the approximate solu- 
tion of the Schrodinger equation and related problems into two main 
categories: 

1. Methods with constant coefficients 

2. Methods with coefficients depending on the frequency of the prob- 
lem ^. 

The purpose of this paper is to introduce a new methodology for the 
construction of numerical methods for the approximate solution of the 
one-dimensional Schrodinger equation and related problems. The new 
methodology is based on the requirement of vanishing the phase-lag and 
its derivatives. The efficiency of the new methodology will be studied 
via the error analysis and the application of the investigated methods 
to the numerical solution of the radial Schrodinger equation. 

More specifically, we will develop a family of hybrid Numcrov-type 
methods of sixth algebraic order. The development of the new family is 
based on the requirement of vanishing the phase-lag and its derivatives. 
We will investigate the stability and the error of the methods of the 
new family. Finally, we will apply both categories of methods the new 
obtained method to the resonance problem. This is one of the most 
difficult problems arising from the radial Schrodinger equation. The 
paper is organized as follows. In Section 2 we present the theory of 
the new methodology. In section 3 we present the development of the 
new family of methods. The error analysis is presented in section 4. 
In section 5 we will investigate the stability properties of the new 

^ When using a functional fitting algorithm for the solution of the radial 
Schrodinger equation, the fitted frequency is equal to: + l)/x^ + V{x) — fe^j 
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developed methods. In Section 6 the numerical results are presented. 
Finally, in Section 7 remarks and conclusions are discussed. 

2. Phase-lag analysis of symmetric multistep methods 

For the numerical solution of the initial value problem 

y" = f{x,y) (3) 

consider a multistep method with m steps which can be used over the 
equally spaced intervals {xjj^Q G [a, 6] and h = — Xi\, i = 

0(l)m- 1. 

If the method is symmetric then Oj = am-i and hi = bm-i, i = 

o(i)LfJ. 

When a symmetric 2fe-step method, that is for i = — A;(l)/c, is applied 
to the scalar test equation 

y" = -cu^y (4) 
a difference equation of the form 

Ak{H) yn+k + ... + Ai{H) vn+i + Aq{H) yn + 

+Ai{H)yn-i + ... + Ak{H)yn-k = ^ (5) 

is obtained, where H = uh, h is the step length and Aq{H), Ai{H), . . ., 
Ai;{H) are polynomials of H. 

The characteristic equation associated with (5) is given by: 

AkiH) A'^ + ... + Ai{H) X + Ao{H) + Ai{H) A'^ + ... (6) 

+AkiH)\-'^ = (7) 

THEOREM 1. [102] The symmetric 2k-step method with character- 
istic equation given by ( 6) has phase-lag order q and phase-lag constant 
c given by 

-n -|_ n(H1+4:\ - 2Afc(if) cos(fciJ)+...+2Aj(JJ) cosQ' g)+...+Ao(JJ) 

cn -r^\n )- 2k-^ Ak{H)+...+2p Aj{H)+...+2Ai{H) y'^' 

The formula proposed from the above theorem gives us a direct 
method to calculate the phase-lag of any symmetric 2k- step method. 
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3. The New Family of Numerov-Type Hybrid Methods - 
Construction of the New Methods 

3.1. First Method of the Family 

We introduce the following family of methods to integrate y" = f{x,y) 

yn = yn- Go (y'n+i -'^yn + y'n-i) 

y-n+l + Cl Vn + yn-1 = [^0 (l/n+l + l/n-l) + blVn] (9) 

The application of the above method to the scalar test equation (4) 
gives the following difference equation: 

Al{H) Vn+l + Aq{H) Vn + Ai{H) Vn-l = 

where H = uoh, h is the step length and ^o(-f^) and Ai{H) are polyno- 
mials of H. 

The characteristic equation associated with (10) is given by: 



Ai{H)\ + Aq{H)+Ai{H)\-^={) (10) 



where 



Ao{H) = ci+H^bi-2H'^ bi ao 

By applying A; = 1 in the formula (8), we have that the phase-lag is 
equal to: 

phi = ^(^) + MH) 



2Ai{H) 

1 2 (1 + 6o + 5^ cos{H) + a + H'^ bi - 2H'^ bi qq 



(11) 



2 1 + H'^bo + H'^bi ao 

We demand that the phase-lag is equal to zero and we consider that: 



bo = ^^b, = ^, Cl = -2 (12) 



Then we find out that: 



_-12cos{H)-cos{H)H^ + 12-5H'^ 
~ 10cos{H)H^-10H^ 
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For small values of \H\ the formulae given by (13) are subject to 
heavy cancellations. In this case the following Taylor series expansions 
should be used: 



° 200 5040 144000 4435200 



99066240000 4790016000 
I 3617 , 43867 

592812380160000 250445794959360000 

+ + (14) 

35213055381504000000 ^ ^ 

The behavior of the coefficients is given in the following Figure 1. 
The local truncation error of the new proposed method is given by: 



LTE 



6048 



REMARK 1. The method developed in this section is the same with 
the obtained by Simos in [116] 

3.2. Second Method of the Family 

Consider the family of methods presented in (9). 

The application of the above method to the scalar test equation (4) 
gives the difference equation (10) and the characteristic equation (10). 

By applying A: = 1 in the formula (8), we have that the phase-lag is 
given by (11). The first derivative of the phase-lag is given by: 

. 1 r4-2Tosin(iJ) + 2iJ6i-8iJ3 6iao 
""'=2 

1 {2TQCos{H) + ci + H'^hi-2H'^hiaQ){2HhQ + AH^hiaQ) 

To = l + H^bo + H^bi ao 
T4 = 2{2Hbo + 4H^bi ao) cos(iJ) (16) 

We demand that the phase-lag and its derivative are equal to zero 
and we consider that: 

^ = 4f. = 5 (17) 
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Behavior of the coefficient a 
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Figure 1. Behavior of the coefficient ao of the new method given by (13) for several 
values of H . 



Then we find out that: 



ao 



-sin(F) + 10 + 2 cos{H) H - 12 sm{H) 
10 sm{H) - 40 cos{H) + 40 
ci = (24 cos (2 H) + 24 - 48 cos{H) + H'^ cos (2 H) 
-9H'^ + 8 cos{H) H"^ -6H^ sm{H) 
-12sm(i?) H)/{6sm{H) H - 24cos(i?) + 24) 



(18) 
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For small values of \H\ the formulae given by (18) are subject to 
heavy cancelations. In this case the following Taylor series expansions 
should be used: 



ao 



1 +-!-H^ + ^l^HU 



509 



200 3780 



5443200 



769824000 



+- 



2833543 



88268019840000 



(19) 



+ 



4912333 



3177648714240000 



H 



10 



+ 



288303913 



3889442026229760000 



165095552521 



46556621053970227200000 



H 



14 



+ - 



15619496804053 



92182109686861049856000000 



+ 



18144 16329600 

308851 537907 
H + 



461894400 



H 



12 



105921623808000 



3813178457088000 



H 



16 



(20) 



The behavior of the coefficients is given in the following Figure 2. 



Figure 2. Behavior of the coefficients of the new method given by (18) for several 
values of H. 



The local truncation error of the new proposed method is given by: 



LTE 



18144 



3y(8)+4a;2y(6) + 



(21) 
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3.3. Third Method of the Family 



Consider the family of methods presented in (9). 

The apphcation of the above method to the scalar test equation (4) 
gives the difference equation (10) and the characteristic equation (10). 

By applying A; = 1 in the formula (8), we have that the phase-lag 
is given by (11). The first derivative of the phase-lag is given by (16). 
The second derivative of the phase-lag can be written as: 



We demand that the phase-lag and its first and second derivative 
are equal to zero and we consider that: 



•• IT3-A T2 sm(H) - 2 Ti cos(H) + 2 61 - 24 61 ao H'^ 

P^^=2 

(2 T2 cos(ij") - 2 Ti sin(ij") + 2 iJ bi - 8 iJHi ap) T2 

T? 

(2 Ti cos(iJ) +ci + H^bi-2H'^ bi uq) Ti'^ 
+ ^rs 

1 (2 Ti cos(ij') + ci + H^bi-2H^ bi ap) (2 bp + 12 hi ap H^) 

2 T? 



Ti = l + H^ bp + H'^bi ap 
T2 = 2Hbp + iH^ bi ap 



T3 = 2{2bo + 12biaoH^)cos{H) (22) 




(23) 



Then we find out that: 
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ao = - (cos{H) + 12 cos{H) H - 12 sm{H) + 3 sm{H) H^^ / 

((cos(iJ)2 + 16cos(i?")2 H + 5cos{H) H"^ sin(iJ) 
+72 cos{H) sui{H) + 2 cos{H) + 32 cos{H) H 
+2 sm{H) H'^ -48H -2H^ -72 sm{H)^H^^ 

ci = ^ (24 cos(ii")2 + cos{Hf + 96 cos{Hf 

+cos{H) sm{H) + 12 cos(i7) iJ^ 
-24 cos(i7) sin(i/) i7 - 96 cos(i/) + cos(i7) H'^ 

-sin(il) iJ^ - 2 i?^ - 60 siii(iJ) - 48 i?^) / 

(cos(ii") + 7sin(ii") + 8 - 8cos(/i")) 

61 = --{cosiHf + 16 cos(i/)2 if + 5 cos(if) if^ sin(if) 
6 V 

+72 cos{H) sm{H) + 2 cos(i7) 
+32 cos(ii) H + 2 sm{H) - 4& H - 2 - 72 sin(ii)) / 

{h (cos{H) H^ + 7 sm{H) H + 8-8 cos(/i)) ) (24) 



For small values of \H\ the formulae given by (24) are subject to 
heavy cancellations. In this case the following Taylor series expansions 
should be used: 
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° 200 2520 907200 1197504000 

I 18427 8 _ 669341 ,0 

980755776000 98075577600000 

13764419 .2 281298850211 ^4 
j3 — ■ — M 



25184162304000000 5747730994317312000000 



161773544323 



103459157897711616000000 



1 o 17 in 43 in 

^ 6048 2721600 57480192 

1515133 25819 



23538138624000 4483454976000 
5 1 „R 11 „« 2353 



6 3024 725760 1437004800 

I 186533 12 , 112457 ,4 

1307674368000 8826801984000 
I 1635421 
1440534083788800 ^ ^ 



The behavior of the coefficients is given in the foUowing Figure 3. 
The local truncation error of the new proposed method is given by: 



3.4. Fourth Method of the Family 



Consider the family of methods presented in (9). 

The application of the above method to the scalar test equation (4) 
gives the difference equation (10) and the characteristic equation (10). 

By applying A; = 1 in the formula (8), we have that the phase-lag is 
given by (11). The first derivative of the phase-lag is given by (16). The 
second derivative of the phase- lag is given by (22) . The third derivative 
of the phase-lag can be written as: 
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1 Tg - (") T8 cosf //) + 2 T", siu(//) - -18 hi no H 
phi = ■ 

2 Ts 

3 (2 T7 cos{H) - 4 Tg sm{H) - 2 T5 cos{H) + 2 61 - 24 61 oq H'^) Tg 
2 T? 

3 (2 Tg cos(ff) - 2 T5 sin(if) + 2 5i - 8 ifHi oq) Tg^ 
+ T? 
_ 3 (2 Tg cos{H) - 2 T5 sin(F) + 2 iJ 61 - 8 iJHi ao) T7 _ 3 Tg Tg^ 
~2 T? T5^ 

3T6TgT7 12T66iaog 

T53 Ts^ 
T5 = 1 + HHo + HHi ao 
Te = 2 T5 cos(ii') + ci+H^bi-2H'^ bi uq 
T7 = 26o + 1261 ao -ff^ 
Tg = 2iJ6o + 4iJ^ 61 ao 
Tg = 48 bi ao H cos(iJ) - 6 T7 sin(i^;B7) 



We demand that the phase-lag and its first, second and third deriva- 
tive are equal to zero and we find out that: 
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ao = - (3 cos{Hf + cos{Hf iJ^ + 2 iJ^ - s) 

/(^(6cos{Hf H + 6sm{H) cos{Hf 

-2 cos{Hf sm{H) + cos{Hf 
+3 cos{Hf H-6 cos{H) sm{H) 
-Acos{H) H'^ sm{H) - 12 cos{H) H + 2 

+3H + 12sm{H) H^'^H^ 

ci = -2 (-12 cos{Hf H + cos{Hf - 21 cos{Hf H 
-12sm(fl') cos{Hf - 4cos(iI)2 H'^ sm{H) + 12cos(iJ) sin(//) 
-8 cos(iJ) if^ sin(iJ) + 24 cos(iJ) /i" + 2 + 9 + 24 sm{H) H"^^ / 

{cos{Hf - 21 cos{Hf H + 8 cos(iJ) sm{H) 
-12cos(iJ)ii"- 12cos(i?) sin(if) +4sin(i7) H"^ 
+33 + 12 sm(i7) + 2 iJ^) 

60 = -2 (3 cos{Hf H + cos(iI)2 + 6 cos{H) sm{H) 
+4 cos(ii") H"^ sm{H) + 6 cos(i?) H 
+2 sin(iJ) i?^ - 9 - 6 sin(/f) + 2H^^/ 

( (cos(if)^ - 21 cos{Hf H + 8 cos{H) sm{H) 
-12 cos{H) H-12 cos{H) sm{H) 
+4 sm{H) H'^ + 33H + 12 sm{H) + 2 H^^ H^^ 

61 = 4 (6cos(i/)3 + 6sm(i/) cos(i/)2 - 2cos(i/)2 sm{H) 
+cos{Hf + 3 cos{Hf H-6 cos{H) sm{H) 
-4 cos(iJ) H^ sm{H) - 12cos{H) H + 2 H^ + 3 H + 12 sm{H) H^^ / 

(^(cos{Hf H^ - 21 cos{Hf H + 8 cos(F) iJ^ siii(iJ) - 12 cos{H) H 

-12cos(if) sm{H) + 4sin(if) i/^ + 33if + 12sin(if) + 2H^^H'^^8) 

For small values of \H\ the formulae given by (28) are subject to 
heavy cancellations. In this case the following Taylor series expansions 
should be used: 
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a 1 I 1 jj2 , 29 4 1433 , 

° 200 1260 504000 1164240000 

63101 ^8 2228861 ,q 
H H 



363242880000 127135008000000 

8804897 240953700049 

H H H 



77806624896000000 20489596600 1 1 264000000 



9699610781879 



_| .^^^.^ 

819583864004505600000000 

= + 6048 + 43200 + 532224 
I 41 _ 601 

5943974400 24141680640 

6 ^ ^ jj4 31 
° 12 1008 181440 



221 619 



10 



13685760 1345344000 



25031 .2 84256583 

H H H 



174356582400 2667655710720000 
1030007057 

"^290289444157440000 " ' 

6 504 90720 
3271 ^8 35293 



47900160 4540536000 
36019 „io 47333617 



+ ^14 



87178291200 1333827855360000 
294008389 .f, 

24562952967168000 ^ ' 

The behavior of the coefficients is given in the foUowing Figure 4. 

The local truncation error of the new proposed method is given by: 



LTE = ^ + Au' + 6a;^ + Au' + y„) (30) 



4. Error Analysis 

We will study the following methods: 
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Figure 4- Behavior of the coefficients of the new method given by (28) for several 
values of H. 

— The First Method of the Family (mentioned as PLl) 

— The Second Method of the Family (mentioned as PL2) 

— The Third Method of the Family (mentioned as PL3) 

— The Fourth Method of the Family (mentioned as PLA) 
The error analysis is based on the following steps: 

— The radial time independent Schrodinger equation is of the form 

y"{x) = f{x)y{x) (31) 

— Based on the paper of Ixaru and Rizea [25], the function f(x) can 
be written in the form: 

f{x)=g{x)+G (32) 
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where g{x) = V{x) — Vc = g, where Vc is the constant approxima- 
tion of the potential and G = = Vc — E. 



(i) 

We express the derivatives yn , « = 2, 3, 4, . . . , which are terms of 
the local truncation error formulae, in terms of the equation (31). 
The expressions are presented as polynomials of G. 



Finally, we substitute the expressions of the derivatives, produced 
in the previous step, into the local truncation error formulae. 



Based on the procedure mentioned above and on the formulae: 

= (V(x)-Fe + G)y(x) 

^ = ^^^^^ + ' ^1 

+(V(x)-F, + G)(^y(x)) 

+4(Av(x))V(x) 
+6(V(x)-Fe + G)(^y(x))(^V(x)) 
+4(U(x)-F, + G)y(x)(^V(x)) 
+(V(x)-F, + G)2(^y(x))... 



we obtain the following expressions: 
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The First Method of the Family 



LTEpLi 



1 



1 



g048g(x)y(-)G^ + (-2^(^g(x))y(x) 

d , , d 

^ ■(— g(3^))y(3^)-T^(3z^Tg(^))(x:y(^)) 



1 

336 



2016 ^da;4 

?(x)(^y(x))(^g(x)) 



1512 

^g(x)y(x)( 



dx 

dx"^ 



^(x:g(^))'y(^)- ^ 



1 



252 ^dx 

d6 



6048 ^d^^^^^^^^^ 



1 



y(xfY(x)) G 
2016°^ ; 

xjgW)(^yW) 



1008 ^dx^ 

3^ ^^^^ ^^^^ ^^^^^ - 2^ ^^^^ 



13 . d 



3024 S^")) - 4 ^("^^ 



1 



y(-) (£2 g(-)) - ^ g(-) y(-) (£ g(-)) 

1 



^2 



126^^^("))(l^("^)^^^(")^ 

d 



6048 



da; 

x)V(x) (33) 
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The Second Method of the Family 



= [^^2 + ik ^T. 

+^g(^)V(x))G^ + (^(^g(.))y(.) 

+ 324 ^^"^^ + 378 t ( ^^"^^ 

^ g(x) y(x) g(x)) + ^ g(x))^ y(x) 

1 1 

+^g(-)^y(-))G + ^(^g(x))y(x) 



+ 



+ 



1008 ^^"^^ ^("^^ + 378 ^^"^ ^^^^ S^^^) 



13 , d 



d3 



+ 



252 



d 



+ 2016 ^("^1^' + 3024 S^"^^ ^^"^ 

+ 5U4 ^Tx ^(^^^ 

d^ 11 d^ 

g(^)) + ^ dx? j{x) g(x)) 

+^g(-)y(-)(^g(-))^ + ^~^^^^ 



+ 126fe^(^))(d^^^^))^d.^ 



l{xfy{x)\{M) 
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The Third Method of the Family 



+ + ^4 

1 



+— g{x) y(a;)) (— g(a;)) 



/ 1 , , , , 1 



504 



d 



1008 ^^"^ + 1512 ^ S^"))' 

1 \ 1 

+^g(x)^(-))G+^(^g(.))y(.) 

loUs^^^^^^^^^y^^)) 
+^g(^)y(^)(^g(^)) 



+ 



5 . d /\\9 / \ 1^ 



^3 



2016 ^da;2 



3024 Mx^(")^^(")fe^(")) 
1 d ci^ 

+— g(^)(x:y(^))(x:3g(^)) 



d 



d^ 



252 °' ' 'dx ' ^dx^ 
11 , ,o , , , 



+ 



216 



(x)y(x)(^g(x))2 + ^ 



(35) 
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The Fourth Method of the Family 



^T^- = [(554 + €^ 

^ y(a^) g(-)) + 4 g(x)) G 



+ 3^gWy(^)(/^2.v^.. ■ 504 
+ 6^ S^")^ + lOM S^J)) 
+ 3^ S(^) ^^^^ + ^ 

+ 2^ s^^)) + 554 ^1 

+ ^g(.x)S-(.x)(A_g(x)) 

+ 2l6 ^^^^ ^^''^ ^^"'^^^ + 60i8 ^^""^l ^^^^ 
We consider two cases in terms of the value of E: 

- The Energy is close to the potential, i.e. G = — £' « 0. So only 
the free terms of the polynomials in G are considered. Thus for 
these values of G, the methods are of comparable accuracy. This 
is because the free terms of the polynomials in G, are the same 
for the cases of the classical method and of the new developed 
methods. 

- G S> or G ^ 0. Then | G | is a large number. So, we have the 
following asymptotic expansions of the equations (33), (34), (35) 
and (36). 

The First Method of the Family 

LTEpLi = h' g(^) y(^) G^ + ...) (37) 



Paper_l_SAV.tex; 15/11/2008; 9:59; p. 19 



20 T.E. Simos, Z.A. Anastassi, D.S. Vlachos 

The Second Method of the Family 

LTEPL2 = g(x)) j{x) + g(x)) J{x)) 

+^g(^)'yW)G^ + ---) 

The Third Method of the Family 

LTEPL3 = g(x)) y(x) G' + ...) (38) 

The Fourth Method of the Family 

^T^-' = iiw4 (|r ^(-» ^(^^' + 736 <^ ^<^» 

+ 3^8 ^f"^* + 5S4 ^(^>) « + ■ ■ 

Prom the above equations we have the following theorem: 

THEOREM 2. For the First Method of the New Family of Methods 

the error increases as the third power of G. For the Second and Third 
Methods of the New Family of Methods the error increases as the second 
power of G. For the Fourth Method of the New Family of Methods the 
error increases as the first power of G. It is easy one to see that the 
coefficient of the second power of G in the case of the second method 
of the New Family of Methods is 1.583333333 times larger than the 
coefficient of the second power ofG in the case of the third method of the 
New Family of Methods. So, for the numerical solution of the time in- 
dependent radial Schrddinger equation the new obtained Fourth Method 
of the New Family of Methods is the most accurate one, especially for 
large values of \ G \=\ Vc — E \. 



5. Stability Analysis 

We apply the new family of methods to the scalar test equation: 

y" = -t^y, (40) 
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where 00. We obtain the following difference equation: 

Ai{H,s)yn+i+ Ao{H,s)yn + Ai{H,s)yn-i = 

where s = th, h is the step length and Aq{H,s) and Ai{H,s) are 
polynomials of s. 

The characteristic equation associated with (41) is given by: 

Ai{H, s) s + Ao{H, s) + Ai{H, s) s'^ = (41) 

where 

Ai{H,s) = l + s'^bo + s^bi ao 
Ao {H, s) = ci + 61 - 2 s"^ 61 ao (42) 

DEFINITION 1. (see [19]) A symmetric four-step method with the 
characteristic equation given by (4.1) is said to have an interval of 
periodicity {0,Wq) if, for all w G {0,Wq), the roots Zi, i = 1,2 satisfy 

zi^2 = e^''^"'\\zi\<l,i = 3,4: (43) 
where 6{t h) is a real function of t h and s = th . 

DEFINITION 2. (see [19]) A method is called P-stable if its interval 
of periodicity is equal to (0, 00). 

THEOREM 3. (see [20]) A symmetric two-step method with the 
characteristic equation given by (4I) is said to have a nonzero interval 
of periodicity (O, Sq) if, for all s G (O, Sq) the following relations are 
hold 

Pi{H,s)>0, P2{H,s)>0, (44) 
where H = ujh, s = th and: 

Pi{H,s) = Ao{H,s) + 2Ai{H,s)>0, 
P2iH,s)=Ao{H,s)-2Ai{H,s)>0, (45) 

DEFINITION 3. A method is called singularly almost P-stable if its 
interval of periodicity is equal to (0, 00) — S^ only when the frequency of 
the phase fitting is the same as the frequency of the scalar test equation, 
i.e. H = s. 

^ where iS is a set of distinct points 
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Based on (42) the stability polynomials (45) for the new developed 
methods take the form: 

Pi{H,s) = ci + v'^bi + 2 + 2v'^ bo, 
P2{H,s) = ci + v"^ hi - iv^ biao - 2 - 2v'^ ho (46) 

In Figures 5,6,7 and 8 we present the s — H planes for the methods 
developed in this paper. A shadowed area denotes the s — H region 
where the method is unstable, while a white area denotes the region 
where the method is stable. In Figure 5 we present the s — H plane for 
the first method of the new family of method developed in this paper 
(paragraph 3.1). In Figure 6 we present the s — H plane for the second 
method of the new family of method developed in this paper (paragraph 
3.2). In Figure 7 we present the s — H plane for the third method of the 
new family of method developed in this paper (paragraph 3.3). Finally, 
in Figure 8 we present the s — H plane for the fourth method of the 
new family of method developed in this paper (paragraph 3.4). 

Stability Regions 
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Figure 5. s — H plane of the first method of the new family of method developed in 
this paper (paragraph 3.1) 
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Figure 6. s — H plane of the second method of the new family of method developed 
in this paper (paragraph 3.2) 



In the case that the frequency of the scalar test equation is equal 
with the frequency of phase fitting, i.e. in the case that H = s, we have 
the following figure for the stability polynomials of the new developed 
methods. A method is P-stable if the s — H plane is not shadowed. 
From the above diagrams it is easy for one to see that the interval of 

periodicity of all the new methods is equal to: (0, vr^). 



REMARK 2. For the solution of the Schrddinger equation the fre- 
quency of the exponential fitting is equal to the frequency of the scalar 
test equation. So, it is necessary to observe the surroundings of the first 
diagonal of the w — H plane. 
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Figure 7. s — H plane of the third method of the new family of method developed 
in this paper (paragraph 3.3) 



6. Numerical results - Conclusion 

In order to illustrate the efficiency of the new methods obtained in 
paragraphs 3.1 - 3.4, we apply them to the radial time independent 
Schrodinger equation. 

In order to apply the new methods to the radial Schrodinger equa- 
tion the value of parameter v is needed. For every problem of the 
one-dimensional Schrodinger equation given by (1) the parameter v 
is given by 

v = ,/W{x)\ = ^\V{x)-E\ (47) 
where V{x) is the potential and E is the energy. 

6.1. Woods-Saxon potential 

We use the well known Woods-Saxon potential given by 

Uq UqZ 



with z = exp [{x — Xq) /a] , uq = —50, a = 0.6, and Xq 



(48) 



7.0. 
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Figure 8. s — H plane of the fourth method of the new family of method developed 
in this paper (paragraph 3.4) 
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Figure 9. Stability polynomials of the new developed methods in the case that H = s 



The behavior of Woods-Saxon potential is shown in the Figure 10. 
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Figure 10. Behavior of the coefficients of the new method given by (24) for several 
values of H . 



It is well known that for some potentials, such as the Woods-Saxon 
potential, the definition of parameter v is not given as a function of x 
but it is based on some critical points which have been defined from 
the investigation of the appropriate potential (see for details [13]). 

For the purpose of obtaining our numerical results it is appropriate 
to choose V as follows (see for details [13]): 



V-50 + £;, 

V -37.5 + E , 
s/-2b + E, 
V-12.5 + -E, 



forxG [0,6.5-2h], 
for X = 6.5 — h 
for X = 6.5 
for X = 6.5 + h 
for X G [6.5 + 2h, 15] 



(49) 
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6.2. Radial Schrodinger Equation - The Resonance 
Problem 

Consider the numerical solution of the radial time independent Schrodinger 
equation (1) in the well-known case of the Woods-Saxon potential (48). 
In order to solve this problem numerically we need to approximate 
the true (infinite) interval of integration by a finite interval. For the 
purpose of our numerical illustration we take the domain of integration 
as X € [0,15]. We consider equation (1) in a rather large domain of 
energies, i.e. E e [1, 1000]. 

In the case of positive energies, E = k"^, the potential dies away 
faster than the term ^^^^ and the Schrodinger equation effectively 
reduces to 

y"ix) + - ii^) yix) = (50) 

for X greater than some value X. 

The above equation has linearly independent solutions kxj[{kx) and 
kxni{kx) where ji{kx) and ni{kx) are the spherical Bessel and Neu- 
mann functions respectively. Thus the solution of equation (1) (when 
X — GO ) has the asymptotic form 



y{x) 2± Akxji{kx) — Bkxni{kx) 
sin { kx — + tanSicos ( kx 



~ AC 

where Si is the phase shift, that is calculated from the formula 



(51) 



_ yix2)S{xi) - yixi)S{x2) , . 

- yix,)C{x,)-y{x2)C{x,) ^^^> 

for xi and X2 distinct points in the asymptotic region (we choose xi 
as the right hand end point of the interval of integration and X2 = 
xi—h) with S{x) = kxjiikx) and C{x) = —kxni{kx). Since the problem 
is treated as an initial-value problem, we need yo before starting a 
one-step method. Prom the initial condition we obtain yo- With these 
starting values we evaluate at xi of the asymptotic region the phase 
shift Si- 

For positive energies we have the so-called resonance problem. This 
problem consists either of finding the phase-shift 5i or finding those E, 
for E e [1, 1000], at which Si = ^. We actually solve the latter problem, 
known as the resonance problem when the positive eigenenergies lie 
under the potential barrier. 

The boundary conditions for this problem are: 

y(0) = 0, y{x) = cos ^V^a;^ for large x. (53) 
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Resonance Problem. E= 163. 215341 


▲ — 


Method 1 


A — A- 


-A Method 1 1 


• • 


-• Method 1 1 1 





-© Method IV 


B— H- 


-□ Method V 


□ □ 


-H Method VI 


♦ — ♦- 


Method VI 1 




Number of Function Evaluations (NFE) 



Figure 11. Error Errmax for several values of n for the eigenvalue Ei — 163.215341. 
The nonexistence of a value of Errmax indicates that for this value of n, Errmax is 
positive 



We compute the approximate positive eigenenergies of the Woods- 
Saxon resonance problem using: 

— The Numerov's method which is indicated as Method I 

— The Exponentially-fitted four-step method developed by Raptis 
[16] which is indicated as Method II 
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Resonance Problem. £=989.701916 
▲ — ^ — ▲ Method I 
A — A — A Method 1 1 
■— ■ — ■ Method III 
H — H — □ Method IV 
4 — — ^ Method V 
^ — e-^ Method VI 
♦ — ♦ — ♦ Method VI I 




Number of Function Evaluations (NFE) 



Figure 12. Error Errmax for several values of n for the eigenvalue E-i — 989.701916. 
The nonexistence of a value of Errmax indicates that for this value of n, Errmax is 
positive 



— The Two-Step Numerov-type Method with minimum phase-lag 
produced by Chawla and Rao [23] which is indicated as Method 
III 



— The new Two-Step Numerov-Type Method with phase-lag equal 
to zero obtained in paragraph 3.1 which is indicated as Method 
IV. 
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— The new Two-Step Numerov-Type Method with phase-lag and its 
first derivative equal to zero obtained in paragraph 3.2 which is 
indicated as Method V. 

— The new Two-Step Numerov-Type Method with phase-lag and its 
first and second derivatives equal to zero obtained in paragraph 
3.3 which is indicated as Method VI. 

— The new Two-Step Numcrov-Typc Method with phase-lag and 
its first, second and third derivatives equal to zero obtained in 
paragraph 3.4 which is indicated as Method VII. 

The computed cigenenergies arc compared with exact ones. In Figure 
11 we present the maximum absolute error logio (Err) where 

Err = lEodiculated ^accuratel (^^) 

of the eigenenergy Ei, for several values of NFE = Number of Func- 
tion Evaluations. In Figure 12 we present the maximum absolute error 
logiQ {Err) where 

Err = \E(,(ii(.y^ig^fg^ -2'accurote I (^^) 

of the eigenenergy E^, for several values of NFE = Number of Function 
Evaluations. 



7. Conclusions 

In the present paper we have developed a family of methods of sixth 
algebraic order for the numerical solution of the radial Schrodinger 
equation. 

More specifically we have developed: 

1. A Two-Step Numerov-Type Method with phase-lag equal to zero 

2. A Two-Step Numerov-Type Method with phase-lag and its first 
derivative equal to zero 

3. A Two-Step Numerov-Type Method with phase-lag and its first 
and second derivatives equal to zero 

4. A Two-Step Numerov-Type Method with phase-lag and its first, 
second and third derivatives equal to zero 
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We have applied the new method to the resonance problem of the 
radial Schrodinger equation. 

Based on the results presented above we have the following conclu- 
sions: 

— The Exponentially-fitted four-step method developed by Raptis 
[16] (denoted as Method II) is more efficient than the Numerov's 
Method (denoted Method I). 

— The Two-Step Numerov-type Method with minimum phase-lag 
produced by Chawla and Rao [23] (Method III) is more efficient 
than the Exponentially-fitted four-step method developed by Rap- 
tis [16] (Method II) for the energy 163.215341 and less efficient for 
the energy 989.701916. 

— The new developed methods are much more efficient than the older 
ones. 

— The Two-Step Numerov-Type Method with phase-lag and its first 
derivative equal to zero (Method V) is more efficient than the 
New Two-Step Numerov-Type Method with phase-lag equal to 
zero (Method IV) 

— The Two-Stcp Numerov-Type Method with phase-lag and its first 
and second derivatives equal to zero (Method VI) is more efiicient 
than the Two-Step Numerov-Type Method with phase-lag and its 
first derivative equal to zero (Method V) 

— The Two-Step Numerov-Type Method with phase-lag and its first, 

second and third derivatives equal to zero (Method VII) is more 
efficient than the Two-Step Numerov-Type Method with phase-lag 
and its first and second derivatives equal to zero (Method VI) 

All computations were carried out on a IBM PC-AT compatible 
80486 using double precision arithmetic with 16 significant digits accu- 
racy (IEEE standard). 
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